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A meta-model (or a surrogate model) is the modern name for what was traditionally 
t-H called a response surface. It is intended to mimic the behaviour of a computational 

model M (e.g. a finite element model in mechanics) while being inexpensive to 
evaluate, in contrast to the original model M. which may take hours or even days 
of computer processing time. In this paper various types of meta-models that have 
been used in the last decade in the context of structural reliability are reviewed. More 
specifically classical polynomial response surfaces, polynomial chaos expansions and 
kriging are addressed. It is shown how the need for error estimates and adaptivity in 
their construction has brought this type of approaches to a high level of efficiency. A 
new technique that solves the problem of the potential biasedness in the estimation 
of a probability of failure through the use of meta-models is finally presented. 

Keywords: Structural reliability, uncertainty quantification, meta-models, surrogate 
models, polynomial chaos expansions, kriging, importance sampling. 
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^h 1. Introduction 

p^i The coupling of mechanical models and probabilistic approaches has gain a lot of 

interest in the literature in the last 30 years. More specifically the field of structural 
reliability has emerged in the mid 70's in order to provide methods for assessing 
structures by accounting for uncertainties in both the models and the parameters de- 
scribing the structures (e.g. geometrical parameters, material properties and applied 
forces). Structural reliability is nowadays a mature field with industrial applications 
in domains ranging from civil, environmental, mechanical & aero-space engineering 



c3 
0\ 



(Ditlevsen and Madsen 1996 1. It has given a sound basis for the semi-probabilistic 
structural codes such as the Eurocodes and has lead to performance-based engineering. 

In the last 10 years the field of stochastic spectral methods has blown up based on 
the pioneering work by |Ghanem and Spanos| ( [l991) . It has become a field in itself 
which is at the interface of computational physics, statistics and applied mathematics 
as shown in the literature on this topic, which is published in scientific journals of 
these various domains, namely Prob. Eng. Mech., Structural Safety, Reliab. Eng. Sys. 
Safety, }. Comput. Phys., Comput. Methods Appl. Mech. Engrg., Comm. Comput. Physics 
on the one hand, SIAM }. Sci. Comput., SIAM }. Num. Anal, }. Royal Stat. Soc. on the 
other hand, among others. 

The large majority of computational methods associated with structural safety and 
uncertainty quantification rely upon repeated calls to the underlying computational 
model of the structure or system. For instance Monte Carlo simulation is based on the 
sampling of the input parameters according to their distribution, and the evaluation 
of the model response (or system performance in the context of reliability) for each 
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realization. The large number of calls that is required for precise predictions is usually 
not compatible with costly computational models such as finite element models, even 
when high-performance computing platforms are at hand. This has opened a new 
field of research broadly called meta-modelling. This is the goal of this paper to provide 
some overview of meta-modelling techniques with a focus on their use in structural 
reliability. 

The paper is organized as follows. Section 2 describes the ingredients of uncertainty 
quantification and recalls the limitations of Monte Carlo simulation for structural 
reliability analysis. Section 3 describes the general philosophy of meta-modelling 
and addresses in details three types of meta-models: polynomial response surfaces, 
polynomial chaos expansions and kriging. Section 4 shows how a meta-model can be 
used in the context of importance sampling in order to provide unbiased estimators 
of a probability of failure. Finally Section 5 gathers some application examples. 

2. Problem statement 

2.1 . Computational model 

Let us consider a mechanical system whose behavior is modelled by a set of governing 
equations, e.g. partial differential equations describing its evolution in time. After 
proper discretization using e.g. a finite element (resp. finite difference) scheme, and 
using some suitable solving scheme, the computational model may be cast as follows: 

y = M(x) (1) 

In this equation x 6 T>x G 1R M is a vector describing the input parameters of the 
model. It usually gathers the parameters describing the geometry of the system, the 
constitutive laws of the materials and the applied loading. Vector y 6 M9 gathers the 
response quantities which may contain: 

• the displacement vector or selected components of the latter; 

• the components of the strain (resp. stress) tensor at specific points; 

• internal variables (e.g. plastic strain, damage variables, etc.); 

• a combination of the latter, at a specific point-in-time or at various time 
instants. 

In the sequel the computational model A4 is considered as a black box, which is only 
known point-by-point: if a given set of input parameters xq is selected, running the 
model provides a unique response vector. It is also assumed that this model is purely 
deterministic: running twice the model using the same input vector will yield exactly 
the same output. 

Note that in practice evaluating such a computational model may be almost instan- 
taneous if some analytical solution to the constitutive equations exists, whereas it can 
take hours on high performance computers if it results from a large size finite element 
model (or a workflow of chained models). 

In the sequel the various methods reviewed for uncertainty quantification and 
reliability analysis consider that the model cannot be modified by the analyst but only 
run for a set of input vectors. These methods are termed non intrusive in the context 
of uncertainty propagation. 
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2.2. Probabilistic model 

Let us consider that the uncertainties in the model input parameters are modelled 
by a random vector X with support T>x £ K. and prescribed probability density 
function fx(x). It is beyond the scope of this paper to describe in details how such a 
probabilistic model can be built from the available information (resp. data). However 
the following guidelines may be found in the literature: 

• When some natural variability of a parameter is evidenced through a set of 
measured values, the classical approach consists in using statistical inference 
methods in order to fit the best distribution selected among one or several 
families (e.g. Gaussian, lognormal, Gamma, Weibull, etc.) using e.g. the 
maximum likelihood principle. The goodness-of-fit shall be checked using 
appropriate tests (jStuart et al. |1999|. Eventually the best distribution may 



be selected using criteria such as the Akaike (resp. Bayesian) information 
criterion ( |Akaike| |1973] |Schwartz| [1978) . 

In contrast when no data is available, expert judgment should be resorted to. 
The available information may be "objectively" modelled using the principle 
of maximum entropy <[Jaynes 1982} Kapur and Kesavan 1992b. Guidelines such 



as the JCSS probabilistic model code (Vrouwenv elder 1997) are available in the 



literature for modelling loads and material properties in civil engineering (see 
updated versions of the code at http : / /www .jess, byg . dtu . dk). 

• In situations where few data is available, a prior expert judgment can be 
combined with measurements using the framework of Bayesian statistics. 

2.3. Structural reliability 

In structural reliability analysis the performance of the system is mathematically 
described by a failure criterion which depends on: 

• the (uncertain) mechanical response of the system, say JA(X); 

• possibly additional deterministic parameters d (e.g. a codified threshold) or 
random variables X' (e.g. uncertain resistance) 

The failure criterion is mathematically represented by a limit state function (also called 
performance function) which is conventionally defined as follows: 

• The set of parameters {x, x', d} such that g (M (x),x' , d) > defines the safe 
domain V s . 

• The set of parameters \x,x' , d} such that g (M. (x), x', d) < defines the failure 
domain T>r. 

• The limit state surface corresponds to the zero-level of the g-function. 

Gathering all the parameters into a single notation x 6 T>x C JR. for the sake of 
simplicity, the probability of failure of the system is defined by: 



P / = P[g(X)<0]= / fx(x)dx (2) 

The main difficulty in evaluating the probability of failure in Eq.Q leads in the fact 
that the integration domain is defined implicitely. Moreover the dimension of the 
integral is equal to the number of uncertain parameters M which is usually large. 
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2.4. Classical computational methods 

2.4.1. Monte Carlo simulation 

Monte Carlo simulation is the basic and universal approach to solving the reliability 
problem (Rubinstein and Kroese 2008). Recasting Eq.d2l as: 



R A 



lv f (x)fx(x)dx = ~E 



!©/(*) 



(3) 



the probability of failure is equal to the expectation of the indicator function of the 
failure domain, which may be given the following estimator using a N-sample X = 

<x( k ', k = 1, . .. ,N> made of N independent copies of X: 



P f-N 



i>,( x(t) ) 
jt=i 



N 



(4) 



where Nf is the number of samples that fall into the failure domain. This estimator is 

P f (l-P f )/N. 



-f 



unbiased and mean-square convergent, since its variance Var 

The slow convergence rate <x 1 / yN makes the approach particularly inefficient. It is 
easily shown that the coefficient of variation of the estimator reads: 



CV Pf 



Var P. 



7 



E 



7 



NP 



(5) 



f 



From the above equation one can see that a typical evaluation of a probability of 
failure of the order of magnitude 10 _r with CVp < 10% requires about W r+2 

simulations. This number is not affordable for low probabilities (10~ 3 — 10~ 6 ) as soon 
as the computational cost of each evaluation of g (which includes a run of Ai) is non 
negligible, e.g. when finite element analysis is involved. 

2.4.2. Beyond Monte Carlo simulation 

Various methods have been proposed in the past 30 years in order to solve the 
reliability problem efficiently. Broadly speaking they can be classified as follows: 

• methods that aim at decreasing the computational cost after introducing some 
approximation in the reliability estimation. The First Order (resp. Second 
Order) reliability methods (FORM/SORM) are well established in this area, 
Hasofer and Lind| <|1974); |Rackwitz and Fiessler| <|1978); |Breitung| <|1984|>; 



see 



Kiureghian and de Stefano 



(1991) among others. 



Hohenbi chler et al.|(|1987y7[B reitung (1989); Der Kiureghi an et al.| ( |1987| );|Der 



methods that are derived from Monte Carlo simulation with the goal of im- 
proving the convergence: directional simulation I Ditlevsen et al... 1987; Bjerager 
1988), importance sampling (jHohenbichler and Rackwitz 1988; Melchers.. 1990 
Maes et al.l|1993 |Au and Beck 1999|) and more recently subset simulation (|Au] 



and Beck. 2001 



Kataf ygiotis and Cheung 2005 Hsu and Ching} |2010[ | and 
line sampling (Pradlwarter et al. 200 7) . 

Methods that rely on the use of a surrogate model M. which is fast to evaluate 
and may be used in place of the original model M. for reliability analysis 
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A detailed presentation of the various classical methods can be found in the textbooks 
by pitlevsen and Madsen| ( [T996) ; |Melchers| ( |1999) ; |Lemaire| ( |2009[ . This is the aim of 



the present paper to address the last point and review different classes of surrogate 
models and their use in structural reliability. 

3. Meta-models for structural reliability 

3.1. Introduction 

As introduced above, a meta-model is an analytical function with the following prop- 
erties : 

• it belongs to a specific class of functions (the "type" of the meta-model) and it 
is fully characterized by a set of parameters once the class is selected; 

• it is fast to evaluate: carrying out some large size Monte Carlo simulation on 
the meta-model will be affordable; 

• it is fitted to the original model (also called "true" model in the sequel, e.g. the 
limit state function g) using a set of "observations" of the true model, i.e. a 
collection of input /output pairs (each observation is a computer experiment 
in the present context): 

X={(x®,y®=g(x®)),i = l,...,N} (6) 

In the sequel, we will review the following classical types of meta-models: 

• linear (resp. quadratic) polynomial response surfaces 

• polynomial chaos expansions 

• Gaussian processes (also known as kriging surrogates) 

Note that support vector machines have been recently introduced in the field of struc- 
tural reliability by Hurtado (|2004a E). Coming from the world of statistical learning 



this technique is well adapted to the classification of a labelled population, e.g. a 
set of "failure" (resp. "safe") points in the context of reliability. It has been recently 
combined with subset simulation to provide a highly efficient way of assessing small 
probabilities of failure (Deheeger and Lemaire 2007; Deheeger, 2008; Bourinet et al. 



|2011) , see also |Basudhar and Missoum| (2008l; Basu dhar et al . (2008). The detailed 



presentation of this approach is beyond the scope of the present paper. 
3.2. FOSM and FORM method viewed as linear response surfaces 



The first-order second moment method (FOSM) (Cornell| 1969 1 may be interpreted as 



a type of linear response surface. Indeed Cornell's reliability index is defined as the 
ratio between the mean value and standard deviation of the safety margin defined by 
the performance function g: 

uJX) 

Pc = , g (7) 

V / Var[g(X)] 

The latter variance is then obtained from a Taylor series expansion of the performance 
function around the mean value of the input vector denoted by fi, which reads: 

g{x) = g(l*) + Vg J {fi)-{x- F ) + o(\\x-}i || 2 ) (8) 
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From the above equation one gets : 

Var \g(X)] « E \[g(X) - g( F )] 2 } = X7g J (ji) ■ C ■ X7g (jt) 



where: 



C = IE [(X - y) (X - ji) 1 



(9) 



(10) 



is the covariance matrix of X. In other words, Cornell's reliability index is implicitly 
based on a linear response surface built up around the input parameters' mean values. 
Due to a well-known lack of invariance when changing the performance function 
(Ditlevsen and Madsen 1996 Chap. 5), the famous Hasofer-Lind reliability index B^i 
has been proposed (Hasofer and Lind [1974), whose derivation may be summarized 
as follows in the context of the so-called first-order reliability method (FORM): 

• the input variables are transformed into a standard normal space by a suitable 
isoprobabilistic transform T : x \- > u = T{x); 

• the reliability index 6 m is defined as the algebraic distance between the 
origin of this space and the (transformed) limit state surface; 

• the associated probability of failure is obtained from the linearization of the 
limit state surface at the design point u*: 



P/JBORM = <S>(-Phl) 



(11) 



where: 

B HL =si S ng(T- 1 (0) 



M *=argmin{|| M || 2 : g(T^(u))<o\ 

(12) 



ueTR A 



In this respect again, the probability of failure is evaluated after constructing a 
particular linear response surface, namely that obtained by a Taylor series expansion 
of the transformed limit state (in the standard normal space) around the design point 



u\ 



3.3. Quadratic response surfaces 

Using a more classical setting, quadratic polynomial response surfaces may be defined 
as follows: 



M M 

g(x) = a + ]P Hi x { + Y, an X? + £ a tj x t Xj 

t=l i=l l<i<j<M 



(13) 



which may be condensed as follows: 

10*0 = /W T ■ « 
a 1 = (a , ai, ... ,a M , a n , ... ,a MM , {«,■,■, 1 <i <j< M}) (14) 

f{x) J = (l, x\, ...,x M ,x\, ...,x 2 M , {xjXj , 1 < i < j < M) 



In order to fit the response surface a set of observations 
X — I [x^',g (xW J J , i = 1, . . . ,N > is selected and the vector of coefficients a (of 
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size P) is computed by minimizing the least-square error between g and g: 



N r 



<*ms = arg mm ^ 
The solution to this problem reads: 



,(0\ _ fT/„(j) 



/'(* w )-« 



a MS = ( F T F) F T r 



where F is the information matrix of size N x P : 



//(* W ) 



g(* (1) ) g(* (N) ) 



(15) 



(16) 



(17) 



The use of quadratic polynomial response surfaces as a surrogate of the limit state 
function has been pioneered by Faravelli ( |1989[ |. Various variants have been proposed 
throughout the 90 s depending on the choice of the polynomials (e.g. taking the 
cross terms XjXj into account or not) and the experimental design used for the model 



fitting, see Bucher and Bourgund"|( |1990) ; Rajashekhar and Ellingwood (1993); Kim and 
Na ( |1997) ; Das and Zheng I pOOO) . The use of these approaches together with finite 
element models has been popularized by Lemaire (1998) and |Pendola et al] ( |2000) 
where the response surfaces are refined in an adaptive manner around the design 
point obtained at each iteration (see also Gayton et al. d2003)). Further applications 
may be found in puprat and Sellier ( 2006} (concrete structures) and Leira et al. ( |2005[ >, 
among others. 

As a conclusion, quadratic response surfaces have been widely used in the last 
20 years for reliability analysis. However they lack of versatility in the sense that a 
second order polynomial function can only mimic models with smooth behaviours. 
Moreover it is implicitely supposed that there is a single design point as in FORM, 
around which the approximation may be built. This condition is seldom encountered 
in industrial problems, which has lead to the use of more versatile meta-models. 



3.4. Polynomial chaos expansions 

3.4.1. Some history 

Polynomial chaos (PC) expansions have been introduced in the literature on stochastic 
mechanics in the early 90's by |Ghanem and Span os (1991) and have been limited to 
solving stochastic finite element problems throughout the 90's. In the original setting, 
a boundary value problem is considered in which some parameters are modelled 
by random fields. The quantities of interest are the resulting stochastic displacement 
and stress fields. Thus the use of PC expansions has been intimately associated 
with spatial variability and considered as a separate topic with respect to structural 
reliability for a while. 

Considering the expansion in itself as a meta-model that is suitable for reliability 
analysis has been originally explored by Sudret and Der Kiureghian ( |2000 2002} ; |Su^| 
dret et al. ( |2003) . Later on, the use of PC expansions has blown up with the emergence 
of so-called non intrusive methods. More specifically the regression approach has been 
developed and applied to reliability analysis in Berveiller et al. (2004); Choi et al. 
( |2004| ; Berveiller et al.| ( |2006) , among others. 
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3.4.2. Polynomial chaos basis 

Without going into too much mathematical details, one may consider polynomial 
chaos expansions as an intrinsic representation of a random variable that is defined 
as a function of the input random vector X . In the context of structural reliability the 
limit state function leads to define the "random margin" G = M. (X). The probability 
of failure is then defined by Pt = P (G < 0). Assuming that this variable has a finite 
variance and that the input parameters in X are independent (for the sake of simplicity 
in this presentationFJ, the following representation holds ( poize and Ghanem} |2004): 

G=g(X)= £ a a Y a (X) (18) 

a£]N M 

In this equation the Y a (X) are multivariate orthonormal polynomials in the input vari- 
ables and a tt are coefficients to be computed. Since the component of X are inde- 
pendent, the joint PDF is the product of the margins. For each marginal distribution 
fx (xi) a functional inner product is defined: 

{<h.'<h)i- I <Pi( x ) <p2(x) fxAx) dx (19) 

For each variable i — 1, . . . , M a family of polynomials is then built which satisfies 
the following orthogonality properties: 

(*f' p k ] ) = j v , p® W p k ] (*) hi (*) dx = «j ¥ ( 20 ) 

where fo is the Kronecker symbol which is equal to 1 if j = k and otherwise. 
The norm of polynomial P- is ah which is usually not equal to 1. Thus in order to 
build an orthonormal family, the above polynomials are rescaled. Classical families of 
polynomials correspond to classical types of PDFs, namely Hermite polynomials are 
orthogonal w.r.t to the Gaussian PDF, Legendre polynomials w.r.t to the uniform PDF, 
etc. (|Xiu and Karniadakis 2002 1, see Table [T] for their expressions and the associated 



normalization. 

Once the univariate orthonormal polynomials are available, the multivariate poly- 
nomials are built by tensorization. To each M-tuple a = {&\, . . . ,ocm} £ 1N M one 
associates the polynomial Y a (x) as follows: 

M 

Y,(*)=nYi?(*o (2i) 

The family of Y a 's naturally inherits from the orthonormality of univariate polynomi- 
als so that: 

E[Y tf Y j8 ]=<J, j8 (22) 

w^here S K a is equal to 1 if the M-tuples a and /3 are identical and zero otherwise. 
Once the basis is built, a troncature scheme has to be selected in order to carry out 



a It is always possible to transform the original vector into independent variables using e.g. the Nataf or 
Rosenblatt transform. 
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Table 1 . Classical orthogonal polynomials 



Distribution 



Uniform 


1 


Gaussian 




Gamma 


X 


Beta 


1 



PDF 



-Wl 



2/r 



W/2 



-x z /2 



-U[ 



(l-x)"(l+x) b 
B(a)B(b) 



Orthogonal 
polynomials 



Legendre 
Hermite 

He k (x) 

Laguerre Lf(x) 
Jacobi/f(x) 



Orthonormal 
basis 

H ek (x)/Vk\ 



/£' {x)/Z a ,b,k 



i 1 - 



2 a+i+i T{k+a+l)T{k+b+l) 
7k+a+b+l T(k+a+b+l)T(k+l) 



the computation of the coefficients. The classical setting consists in selecting all the 



M 



polynomials of total degree |«| = ),&{ not greater than a given p, i.e. 

;=i 



g(X)*g^(X) 



oteA 



where A 



{' 



eN 



M 



a < 



(23) 



This type of a priori truncation is somehow arbitrary although a value of p = 2 usually 
provides fair results for estimating the mean and variance of the margin G whereas 
p = 3 is required in order to compute probabilities of failure downto 10 -4 with a 
satisfactory accuracy (Sudret. 20071 . Note that error estimates have been recently 
proposed together with adaptive algorithms in order to avoid the problem of the a 
priori selection of A, as shown in Blatman and Sudret (2010a, 2011 1 

3.4.3. Computation of the coefficients 

The original approach to computing the coefficients of a truncated PC expansion in 
computational stochastic mechanics is of Galerkin-type (Ghanem and Spanos |1991) 
and it is termed intrusive since it requires the ad-hoc derivation of a weak formulation 
of the underlying mechanical problem and its discretization. Non intrusive methods 
have emerged in the literature as of 2002 and may classified as follows: 

• projection methods that make a direct use of the orthogonality properties of 
the PC basis: 



««=E[GT 4 (X)] 



T> x 



g(x)Y K (x)f x (x)dx 



(24) 



The latter integral may be computed using either tensorized or sparse quadra- 
ture rules (G hiocel and Ghanem} |2002 Le Maitre et al. 2002 Keese and 
Matthies, 20051 



stochastic collocation methods, which are based on Lagrange interpolation 
in the stochastic space and are essentially equivalent to the former ( |Xiu and] 
Hesthavenj [2005) |Xlu]|2009) ; 
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regression methods which were introduced in the field of PC expansions in 
(2006) and recently improved through error estimation and 



Berveiller et al. 



adaptivity by B 



atman 



( p009) ; patman and Sudret| ( |2010a}[20lT) . 



In this paper we will concentrate on the latter approach which is similar in essence 
to building a quadratic polynomial response surface as shown in Section 3.3 except 
that the basis functions are now the members of the truncated PC expansion. The 
key idea consists in considering the random margin G as the sum of a truncated PC 
expansion and a residual. 

G=g(X)= £«„Y a (X)+e 



(25) 



The coefficients a a are obtained by minimizing the mean square residual: 

\ 21 



arg min E 

flGR card.4 



g(X) 



£> a T a (X) 
ueA 



which is approximated by using an experimental design of size N as in Eq.|[15): 

-\2 



N 

«MS = arg min V 

a£]R cardyl . =1 



(*«)_ 2>.Y,(*«) 
ueA 



(26) 



(27) 



This leads to solving a linear system as in Eq.|[T6|. In contrast to the projection and 
the stochastic collocation methods which make use of sparse grids, the experimental 
design is selected here so as to be space-filling: Latin Hypercube sampling ( |McKay| 
et al. 19791 or quasi-random numbers are used for this purpose. Empirically the size 
of the experimental design is selected as follows: N=2-3 Card^l (jBlatman 2009). 



3.4.4. Adaptive PC expansions and application to structural reliability 

Once the basis is built and the coefficients have been computed the PC expansion is 
treated as a global polynomial response surface and substituted for the "true" limit 
state function for computing the probability of failure: 



! 7 



jPC 



p (g pc < 



fx(x)dx 



(28) 



f V~ ) J{x:gK(.x)<0}- 

As g FC (x) = V a a Y a (X) is polynomial and straightforward to evaluate, crude 

neA 
Monte Carlo simulation may be used to compute P? c , although more advanced 
simulation techniques such as subset simulation may be used. 

The main difference between this approach and the traditional quadratic response 
surfaces are: 

• an arbitrary high degree of polynomials may be used, which allows one to fit 
complex limit state functions beyond a quadratic approximation. 

• error estimates based on cross-validation techniques ( |Stone 1974} are avail- 
able ( |Blatman and Sudret 2010a) which may be coupled with adaptive algo- 
rithms that automatically detect the best sparse PC representation ( Blatman] 
and Sudret 2011) and increase the maximal degree of the PC expansion as 



long as the prescribed admissible error is not attained. 
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in order to avoid overfitting when computing the coefficients, the size of the 
experimental design N is increased automatically so as to always satisfy the 
condition N < 2 — 3 Card A. 

by construction the PC expansion readily provides useful information on the 
moments of the performance G PC = g (X) as well as useful information on 



sensitivity analysis ( Sudret 2008 ) 



A synthesis of such an approach with applications in structural reliability may be 
found in |Sudretet"aI] ( |2011| |. 



3.5. Kriging meta-models 

3.5.1. Introduction 

Building surrogate models in order to reduce the computational cost associated with 
reliability analysis, stochastic or optimization problems has been given the generic 
name of computer experiments in the statistical literature. In this respect kriging has 
emerged in the last two decades as a powerful tool for building meta-models. It has 
not been used in structural reliability until recently though. 

Historically kriging was named after the South African engineer D. Krige who 



initiated a statistical method for evaluating the mineral ressources and reser ves (jKrige 
1951 1. This opened the field of geostatistics later formalized by Matheron (1963J, see 



Cressie| (1993}; Chiles and Delfinerj (|1999). The term kriging has been coined 



also 

in order to honor the seminal work of D. Krige. The basic idea is to model some 
function known only at a finite number of sampling points as the realization of a 
Gaussian random field. In this setting the sampling space is a "physical" two- or 
three-dimensional space. 

Later, Sacks et al. (1989) introduced the key idea that kriging may also be used in 
the analysis of computer experiments in which: 

• the data is not measured but results from evaluating a computer code, i.e. a 
simulator such as a finite element code; 

• the points where data is collected are not physical coordinates in a 2D or 3D 
space, but parameters in an abstract space of arbitrary size M. 

In contrast to polynomial chaos expansions kriging provides a meta-model that does 
not depend on the probabilistic model for the input random vector X. 

3.5.2. Mathematical setting 

The modern setting of kriging for computer experiments (also called Gaussian process 
modelling (Santner et al. 20031) reads as follows. The function-to-surrogate (e.g. 



the limit state function g in the context of structural reliability) is supposed to be a 
realization of a Gaussian process denoted by Y(x, to) defined as follows: 

Y(x,to)=f (x) T a + Z(x, to) (29) 

In this equation / (x) a is the mean of the process, which is represented by a 
set of basis functions {/;, i = 1, ... ,P} (e.g. polynomial functions) and Z (x, to) 
is a stationary zero mean Gaussian process with variance <7y and autocorrelation 
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function^] 



C YY (x, x') 



a Y R(x- 



x',G), 



[X, X 



£V x xV, 



(30) 



In the above equation 6 gathers all the parameters defining C YY . In practice, square 
exponential models are generally postulated: 



M 



R(x- x', 6) = exp I ]£ ■ 




(31) 



although other types of autocorrelation models such as generalized exponentials or 
the Matern kernel may be used (Santner et al. 2003} . 

In order to establish the kriging surrogate a set of computer experiments is run and 

gathered in a vector Y = (g{x^>), . . . ,g(x^ ') ) . The kriging estimator at a given 

point x 6 T>x is by definition a Gaussian random variate Y (x) ~ J\f (jlo (x) , o~o (x)) 
obtained by requiring that it is the best linear unbiased estimator (BLUE) of g(x) 
conditioned to the observations gathered in X. In other words it is obtained as a linear 
combination of the observations and it is unbiased with minimum variance. After 



some rather lengthy algebra the kriging estimator reads as follows (see Santner et al. 
pO^ / [DubmIrg1 ( |20TT[ Chap. 1)): 



u Y (x) = f (x) T a + r {x) T R 1 (T - F U) 
In this equation the following notation r, R et F is used: 

F, 7 



r(x-x®,Q\, i = 1, ... ,N 

R (x (;) - x®, fl) , i = 1, . . . , N, j = 1, 

Mx®), i = l,...,p,j = l N 



,N 



(32) 

(33) 
(34) 
(35) 



On top of the mean prediction given in Eq.d32}, the kriging approach yields the so-called 
kriging variance Co (x) which corresponds to an epistemic uncertainty of prediction 
that is related to the finite size of the available observation data gathered in X. This 
variance reads: 



4(x) 



f(x) J r{x) T 



OF 1 
F R 



i-l 



r (x) 



(36) 



So far the "regression" part of the kriging estimator fix) a in Eq.(32 and the 
parameters 9 in Eq.lpTjl have not been solved for. The classical approach consists in 
building an empirical BLUE by using a likelihood function from the joint (Gaussian) 
distribution of the observations. Indeed by the underlying assumption of kriging, the 

values gathered in Y form a single realization of a Gaussian vector < Y^\ . . . , Y^ N > V. 

Then the likelihood of the observations in Y is maximized with respect to {a, cr Y , 8}- 



b For the sake of clarity the notation w that recalls the randomness of the various quantities is abandoned in 
the sequel. 
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It may be shown that this optimization yields an analytical expression for { a, o~y } as 
a function of d: 

a = (f t R -1 f) F T R- 1 r (37) 

^ = i(r-F«) T R" 1 (r-Ffl) (38) 



In these expressions the dependence R(0) (see Eq.ipljl) has been omitted for the 
sake of clarity. The best fit values of d are eventually obtained from a numerical 
optimization, see e.g. Marrel et al. (2008); Dubourg (2011) for details. 



The great features of kriging compared to the polynomial response surfaces pre- 
sented in Sections I3.3H3.4I are summarized below: 



• The mean kriging estimator given in Eq. | |32) is interpolating the data, mean- 
ing that jio ( a?W ) = g ( x^ 1 ' J and the kriging variance is zero in these points: 

4(*<o)=o 

• The kriging variance cr~{x) is interpreted as a measure of the epistemic 
uncertainty of prediction in each point x. It shall not to be confused with 
the aleatoric uncertainty represented through random vector X that is related 
to the probability measure P [dx] = fx( x ) dx. Thus o~{x) may be used as an 
indicator for adaptively enrich the experimental design and refine the meta- 
model. 

3.5.3. Adaptive kriging and applications to structural reliability 

Kriging has been first used for structural reliability problems in the contributions by 



Romero et al. (2004) and Kaymaz (|2005). In these early papers the experimental design 



is either fixed or enriched "passively", i.e. not taking into account the information 
brought by the kriging variance. 

Enriching sequentially the experimental design by using a criterion that leads to 
adding points in the vicinity of the limit state function has been proposed by|Bichon| 



et al. (2008) under the acronym EGRA (efficient global reliability analysis). The 



authors define an expected feasibility function EFF(x) that provides an indication on 
the vicinity of the current point to the true limit state function. Starting from the 
premise that only points close to the limit state function bring additional information 
to build the surrogate, they maximize this criterion at each step to get the next point 
to be added to the experimental design for the next iteration. 

Recently a similar approach called AK-MCS (for "active kriging + Monte Carlo 
simulation") has been devised by Echard et al. (2011), who propose an epistemic error 



function which is directly based on the kriging variance. The so-called LT-function is 
defined by: 

uw = W|l 

A small value of U(x) means that either the kriging-limit state function (kriging-LSF) 
fly (x) is close to zero (vicinity of x to the surrogate limit state function ) or Co (x) is 
large (large uncertainty prediction) or both. Once a first kriging prediction has been 
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carried out using a small size experimental design, this U-function is evaluated on 
a large-size Monte Carlo sample. Note that it does not require new calls to the g- 
f unction. The point leading to the smallest value of U(x) is added to the experimental 
design for improving the kriging predictor at the next iteration. 

In both approaches, the active enrichment of the experimental design is carried 
out sequentially, i.e. one single point is added from one iteration to the other. This 
may be considered as a weakness since usually the criteria to optimize show several 
local extrema (in other words, there are several candidate points that could be equally 
added at each iteration). Moreover it is nowadays common to have distributed 
computing facilities which enables parallel evaluations of the limit state function for 
a set of values of x. 

Starting from this observation Dubourg (|2011};|Dubourg et al. (|2011[| have proposed 



to define the following probabilistic classification function: 



n{x,t)=V[nx)<t] = <&(' * /J}? ) (40) 



Cy [X] 



In this equation again, V [•] denotes the Gaussian probability measure associated with 
the epistemic uncertainty of kriging and not the aleatoric uncertainty in X (probability 
measure denoted by P [•]). The vicinity of any point x to the kriging-LSF defined by 
{£ : f*y(£) = 0} * s measured by n{x, 0), which shows some similarity with the U- 
f unction in Eq.<[39). From this function a margin of uncertainty 9JI is defined which 
corresponds to "confidence intervals" around the kriging-LSF: 

971 = {x : -ka ? (x) < u ? (at) < +k a ? (x) } (41) 

where A: is a "number of standard deviations", e.g. k = 1.96 for a 95% confidence 
interval. The enrichement criterion is then defined as the probability of being in the 
margin of uncertainty at point x which turns out to be: 

C{x)=P \?(x) e 0tt| = n{x,ka ? (x)) - n(x,-ka ? (x)) (42) 

This quantity C(x) is then considered as a probability density function (up to a 
constant) from which it is simulated using a Markov chain Monte Carlo algorithm 
such as the slice sampling (Neal, 2003). Instead of picking up a single point as in AK- 



MCS or EGRA, this approach yields a large set of points which concentrate in the 
margin (see Figure [TV 

A limited number K of new experimental points (to be added to X) is then obtained 
by clustering the large sample set and added to the experimental design. Convergence 
criteria are established based on the bounds on Pr obtained by the boundaries of the 
margin 9Jt, see Dubourg et al. ( 2011} for details. 



3.6. Conclusion 

Quadratic response surfaces, polynomial chaos expansions and kriging techniques 
have been reviewed as meta-modelling techniques for structural reliability analysis. 
Whatever their respective advantages and drawbacks, all these meta-models are 
usually used as substitutes for the original limit state function, meaning that the proba- 
bility of failure is simply evaluated as in Eq.([28). In principle there is no guarantee that 
the probability of failure evaluated on the surrogate is equal or sufficiently close to the 
that obtained from the surrogate. A way of addressing this issue is now proposed. 
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Fig. 1. Illustration of the adaptive kriging technique (after Dubourg et al.|j2011} ). The original 
sample set is space-filling. In the later iterations, the sample points fall mainly within the margin 
of uncertainty defined by the kriging variance. Clustering is used in each step in order to select 
a few points to enrich the experimental design. 



4. Meta-models as a means for importance sampling 

4.1. Introduction 

As shown in the previous sections meta-models are usually built in a first step of the 
reliability analysis and then substituted for the "true" limit state function. Even if 
some adaptive scheme allows one to ensure the closeness of the original function and 
its surrogate, the unbiasedness of the estimator of Pt based on the surrogate cannot 
be guaranteed. 

In order to solve this problem, |Dubourg et al.| ( |2011| propose to combine the 
well-known importance sampling technique with kriging in order to get unbiased 
estimates of the probability of failure. As shown below, the surrogate is used to build 
a "smart" importance sampling density instead of being substituted as in Eq.(p8). 



4.2. Reminder on importance sampling 

The reason why crude Monte Carlo simulation is inefficient for evaluating probabil- 
ities of failure is that a typical sample will contain a majority of points located in 
the central part of the input distributions whereas the realizations that lead to failure 
are in the tails. The key idea of importance sampling is to use an instrumental density 
which allows one to concentrate the drawn samples in the region of interest. Consider 
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a non zero probability density function h(x) defined on R . Eq.O may be recast as: 



E;, 



lv f {x) T(xy 



(43) 



where the subscript in E/, [•] means that the expectation is taken with respect to 
the probability measure associated with h. The classical importance sampling (IS) 
estimator is built from a sample set drawn from the instrumental distribution, say 

X h = {4\i = l n}: 



The art of importance sampling consists in using an instrumental density that mini- 
mizes the variance of the latter estimator. Rubinstein and Kroese ( 2008 1 shows that 
the optimal density (which actually reduces the variance of the estimator to 0) reads: 

l v (X)f x {x) lv f (x)fx(x) 

h (X) = ~ r , : z , : = T —= (45) 

However this density cannot be used in practice since it depends on the unknown 
quantity of interest Pf. The key idea of meta-model-based importance sampling (called 
meta-IS for short) is to surrogate this optimal instrumental density by kriging. 

4.3. Sub-optimal instrumental density 

Assuming that some kriging meta-model Y of the limit state function is available, 
the indicator function of the failure domain may be replaced by the probabilistic 
classification function (see Eq.d40b) evaluated for t — 0: 

Indeed if the surrogate is accurate in a given point xq then Co(xo) is close to zero 
and pif(xo) is close to g(xo). If the latter is negative then ^y(xo) / '(t^{xq) —> — °° and 
7r(xo) ~ 0. In contrast, if g(xo) > then u^(xo) / o~<f(xo) — > +oo and rz(xo) ~ 1. Thus 
n(x) is a kind of "smoothed" version of the indicator function lp . This leads to the 
meta-IS instrumental density (to be compared with Eq.d45l) : 



~ h( , _ n{x)f x {x) ^(-v ? (x)/o- ? (x))f x (x) 

h \ x ) = 5 = 5 (47) 



where the normalization constant Pt e reads: 



Pfe = I * (-F? (x)/<r ? (x)) fx(x) dx (48) 

JV X 

As an example the optimal instrumental density and that obtained from the above 
procedure is shown in Figure |2] where the quadratic limit state is taken from Der 



Kiureghian and Dakessian (1998). 
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Fig. 2. Comparison of the optimal (left) and the meta-IS (right) instrumental probability density 
functions (after jDubourg et al.|J201 1| ) 



4.4. Meta-IS estimator of the probability of failure 

Substituting for Eq.fiTl in Eq.d43l and using Eq.d48l one gets: 



where: 



I f — M-corr I fs 



Mi 



n(X) 



(49) 



(50) 



is a correction factor that quantifies the error that has been made by substituting the 
probabilistic classification function tc(x) to the indicator function lj> (X). The Monte 
Carlo estimate of the latter then reads: 



P fmetaIS ~ ° icorr P f £ 



where: 



I N corr lj, (*(*)) 
Ncorr Pi 7T(* W ) 



(51) 



N? 



'> 



M- 



£>(*«) 



1=1 



In the above equation, the first sum is computed using N CO rr realizations following the 
meta-IS PDF h w^hile the second sum is computed using N £ realizations of X. 

The estimator in Eq.|51 1 is proven to be unbiased and its variance may be derived 
straightforwardly Dubourg et al. (2011). From a computational point of view, it is 
important to note that Pr E is evaluated at low cost since it only requires the evaluation 
of the probabilistic classification function Eq.(46 1 which is analytical when the kriging 
surrogate is available. Thus N e = 10 5-6 is affordable, leading to an almost exact 
value of Pf £ . In contrast, computing a corr requires the evaluation the "true" indicator 
function of the failure domain. Thus N CO rr shall be limited to a few hundredw in 
practice. 
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4.5. Meta-IS estimator of the probability of failure 

As a conclusion, using meta-models as a tool for defining a quasi-optimal instru- 
mental density for importance sampling solves the problem of possible biasedness 
resulting from a direct substitution. This is a great advantage of meta-IS compared to 
the approaches by Bichon et al. (2008); Echard et al. (|2011[|. 



In order to optimize the efficiency a trade-off between a) the number of calls to 
the g-function (say, N) required for building the meta-model Y and b) the number 
of samples ~N corr shall be found, since the "total cost" is N + N CO rr- A large N will 
provide a very accurate meta-model that will lead to DC C orr ~ 1 at a high computational 
cost though. In contrast a "medium accurate" surrogate will be compensated for by 
a correction factor a corr ^ 1 . In order to balance these two aspects a procedure based 
on cross-validation is proposed by Dubourg et al. (20111 for stopping the iterative 



improvement of the kriging surrogate (Section 3.5.3 1 when Pe £ is sufficiently close to 



Pf- 



5. Application examples 

5.1 . Frame structure - Polynomial chaos expansion (Blat man and Sudret 



2010a) 



Let us consider a 3-span, 5-storey frame structure as the one sketched in Figure [3] Of 
interest is the top floor displacement when the structure is submitted to lateral loads. 
The associated serviceability limit state function reads: 

g(X) = u max - FEMODEL(X) (52) 

where u max is a given threshold, FEMODEL(») is the finite element model (considered 
as a black box) that evaluates the top floor horizontal displacement, and X gathers 
21 correlated random variables which describe the uncertainty in the member's cross 



section, inertia and applied loads, see Blatman and Sudretj ( 2010a I for the complete 
description. The probability of failure is investigated for different values of the 
threshold, namely Umax = 4 — 9 cm. Reference values are obtained by FORM followed 
by importance sampling (500, 000 model evaluations are used to get a coefficient of 
variation less than 1.0% on PA. 

Estimates of the reliability index are computed by post-processing a full third-order 
PC expansion as well as a sparse PC approximation. The latter is built up by setting 
the target approximation error equal to 10 -3 (this is an overall mean square error 
that drives the convergence of the adaptive PC expansion, see the original paper for 
details). The estimates of the various generalized reliability indices f> = —^>~ 1 (Pf) 
are reported in Table [2] 

From the results in Table [2] it is seen that the sparse polynomial chaos expansion 
allows one to evaluate accurately reliability indices up to /3 = 4 with less than 5% 
error. The parametric analysis w.r.t the threshold is obtained for free since the PC 
expansion is computed once and for all. As a whole the analysis has required 450 
calls to the finite element model. Note that it also yields the statistical moments of the 
response as well as sensitivity results at the same cost (see Blatman and Sudret d2010b ) 
for details). 

It should be emphasized that PC expansions provide an approximation whose 
accuracy is controlled globally, i.e. in a least-square sense. It does not guarantee a 
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Fig. 3. Sketch of a 3-span, 5-story frame structure subjected to lateral loads 



Table 2. Frame structure - Estimates of the generalized reliability index 
/3 = -O -1 (Pf) for various values of the threshold displacement 



Threshold (cm) 


Reference 


Full PCE 


Sparse PCE 




pREF 


is 


e (%) 





e (%) 


4 


2.27 


2.26 


0.4 


2.29 


0.9 


5 


2.96 


3.00 


1.4 


3.01 


1.7 


6 


3.51 


3.60 


2.6 


3.61 


2.8 


7 


3.96 


4.12 


4.0 


4.11 


3.8 


8 


4.33 


4.58 


5.8 


4.56 


5.3 


9 


4.65 


4.99 


7.3 


4.94 


6.2 


Relative error 




1 


10~ 3 


1 


10~ 3 


# terms in PCE 




2,024 




138 


# FE runs 


53,240 


3,724 




450 



perfect control of the accuracy in the tail of the distribution of the limit state function 
though. This means that this approach should not be used for very small probabilities. 
Results reported in the literature show that it is robust up to Pf = 10~ 4 ' ~ 5 . 
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5.2. System reliability - meta-IS ( Dubourg\\201iy 



As an illustration consider the following system limit state function ( Waarts 2000) 



g(x) = min 



/3 + (xi - x 2 ) 2 /10 - Oi + x 2 )/V2 
3 + (xi - x 2 ) 2 /10 + [x 1 + x 2 )/V2 

X 1 -X2 + 7/V2 

\ x 1 -x 1 +7/^/2 ) 



(53) 



where the two components of the input random vector X are independent standard 
normal variables. 

Table 3. Meta-importance sampling - system reliability analysis (after |Dubourg| 
<20TT} ) 





Monte Carlo 


Subset simulation 


Meta-IS 


Computational cost 
C.o.V 


172,000 
2.26 10~ 3 

<5% 


284,195 

2.28 10~ 3 

<3% 


40 + 200 
2.38 10~ 3 

<5% 



The reference results are obtained by crude Monte Carlo simulation so as to obtain 
a coefficient of variation less than 5%. They are reported in Table [3] together with the 
results obtained by meta-importance sampling. It is shown that only 40 runs of the 
limit state function are required in order to build a sufficiently accurate surrogate. 
Then only 200 runs of the true model are required using the sub-optimal instrumental 
density in order to compute the probability of failure within 5% accuracy. The 
computational cost is thus decreased by two order of magnitude. 

6. Conclusions 

Meta-modelling techniques have become inescapable in modern engineering since 
they allow one to address real-world reliability problems at an affordable compu- 
tational cost. In this paper, well-established reliability methods such as FOSM and 
FORM/SORM are reinterpreted as basic meta-models. More advanced techniques 
such as polynomial chaos expansions and kriging have been reviewed and show 
promising features in terms of performance and accuracy. 

This paper also emphasizes the need for engineers to catch up with the latest 
technologies developed in computer science and applied mathematics (including 
statistics) in order to further innovate in their respective fields. 
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